__kernel void computeGrainID(__global int * spaceGrainID,
                             __global int * spaceGrainIDPrv,
                             __global float * spaceGrainRo,
                             __global float * spaceGrainRoPrv,
                             __global int * spaceGrainRecryst,
                             __global int * spaceGrainRecrystPrv,
                             int width,
                             int height,
                             float deltaRo,
                             float RoCritical){

    int gIdx = get_global_id(0);
    int gIdy = get_global_id(1);
    int minusX, minusY, plusX, plusY;   //should be uint

    /*Compute periodic index*/
    if(gIdx - 1 < 0) minusX = width - 1;
    else minusX = gIdx - 1;
    if(gIdy - 1 < 0) minusY = height - 1;
    else minusY = gIdy - 1;
    if(gIdx + 1 > width - 1) plusX = 0;
    else plusX = gIdx + 1;
    if(gIdy + 1 > height - 1) plusY = 0;
    else plusY = gIdy + 1;
    /*ENDCompute periodic index*/

    /*Verify if grain is on board*/
    int isOnBoardNeumann = 0;  //0 if false, other number if true
    int activeGrainID = spaceGrainIDPrv[gIdx * width + gIdy];

    /*Vertical*/
    if (spaceGrainIDPrv[gIdx * width + MinusY] != activeGrainID && spaceGrainRecrystRrv[gIdx * width + MinusY] == 0) isOnBoardNeumann = 1; //up
    else if (spaceGrainIDPrv[gIdx * width + PlusY] != activeGrainID && spaceGrainRecrystRrv[gIdx * width + PlusY] == 0) isOnBoardNeumann = 1;  //down
    /*Horizontal*/
    else if (spaceGrainIDPrv[MinusX * width + gIdy] != activeGrainID && spaceGrainRecrystRrv[MinusX * width + gIdy] == 0) isOnBoardNeumann = 1; //left
    else if (spaceGrainIDPrv[PlusX * width + gIdy] != activeGrainID && spaceGrainRecrystRrv[PlusX * width + gIdy] == 0) isOnBoardNeumann = 1;  //right
    /*END Verify if grain is on board*/

//    for (int i = 0; i < width; i ++){
//            for (int j = 0; j < height; j ++){

                /*Update deployment*/
                //this->updateDeplyment(i, j, deltaRo);
                int rand = std::rand() % 4;     //SOLVE THIS rand() and switch int to uint

                if (isOnBoardNeumann == 1){
                    float delta = (1.0 - (float)rand * 0.1) * deltaRo;
                    spaceGrainRo[gIdx * width + gIdy] += delta;
                    spaceGrainRoPrv[gIdx * width + gIdy] += delta;
                }
                else{
                    float delta = (float)rand * 0.1 * deltaRo;
                    spaceGrainRo[gIdx * width + gIdy] += delta;
                    spaceGrainRoPrv[gIdx * width + gIdy] += delta;
                }
                /*END Update deployment*/

                //space[i][j].recInPrvStep = false;

                /*IF grain is on board and ro > ro_critical*/
                /*new embryo*/
                if (spaceGrainRoPrv[gIdx * width + gIdy] > RoCritical && isOnBoardNeumann == 1 && spaceGrainRecrystRrv[gIdx * width + gIdy] == 0){
                    spaceGrainRecryst[gIdx * width + gIdy] = true;
                    spaceGrainID[gIdx * width + gIdy] = ++this->numGrains; //new embryo grain   //THINK HOW SOLVE IT
                    spaceGrainRo[gIdx * width + gIdy] = 0;
                    //space[i][j].recInPrvStep = true;
                    //this->progress++;
                }

                /*IF neighbor was recrystalized in previous step (t-1) and grain[i, j].ro > ro_in_all_ngh*/
                else if (spaceGrainRecrystRrv[gIdx * width + gIdy] == 0 &&
                        spaceGrainRoPrv[gIdx * width + MinusY] < spaceGrainRoPrv[gIdx * width + gIdy] &&    //up
                        spaceGrainRoPrv[gIdx * width + PlusY] < spaceGrainRoPrv[gIdx * width + gIdy] &&   //down
                        spaceGrainRoPrv[MinusX * width + gIdy] < spaceGrainRoPrv[gIdx * width + gIdy] &&  //left
                        spaceGrainRoPrv[PlusX * width + gIdy] < spaceGrainRoPrv[gIdx * width + gIdy])  //right
                {
                    /*Verify if neighmbour grain is recrystalized and return grain id if yes or 0 if not*/
                    //int newNgh = this->recrystalizeGrain(i, j);
                    int newNgh = 0;

                     int nghID[4];
                    int numRecrystGrain = 0;    //number of recrystalized grain neighbour
                    /*Vertical*/
                    if(spaceGrainRecrystRrv[gIdx * width + MinusY] != 0){
                        nghID[numRecrystGrain] = spaceGrainID[gIdx * width + MinusY];  //up
                        numRecrystGrain++;
                    }
                    else if(spaceGrainRecrystRrv[gIdx * width + PlusY] == 0) {
                        nghID[numRecrystGrain] = spaceGrainID[gIdx * width + PlusY];  //down
                        numRecrystGrain++;
                    }
                    /*Horizontal*/
                    else if(spaceGrainRecrystRrv[MinusX * width + gIdy] == 0) {
                        nghID[numRecrystGrain] = spaceGrainID[MinusX * width + gIdy];  //left
                        numRecrystGrain++;
                    }
                    else if(spaceGrainRecrystRrv[PlusX * width + gIdy] == 0) {
                        nghID[numRecrystGrain] = spaceGrainID[PlusX * width + gIdy];   //right
                        numRecrystGrain++;
                    }

                    if (numRecrystGrain == 0) newNgh = 0;
                    else if (numRecrystGrain == 1) newNgh = nghID[0];
                    else newNgh = nghID[rand() % numRecrystGrain];      //SOLVE THIS rand()

                    if(newNgh != 0){
                        spaceGrainID[gIdx * width + gIdy] = newNgh;
                        spaceGrainRecryst[gIdx * width + gIdy] = 1;
                        spaceGrainRo[gIdx * width + gIdy] = 0;
                        //space[i][j].recInPrvStep = true;
                        //this->progress++;
                    }
                }
//            }
//        }
}